<<<<<<< HEAD ======= >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e Tennessee COVID-19 Update <<<<<<< HEAD ======= >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 
library(tidyverse)
# library(rvest)
# library(albersusa)
library(sf)
library(leaflet)
library(leafpop)
library(formattable)
library(ggiraph)

# library(covdata)

# Default ggplot theme
theme_set(theme_bw() +
            theme(axis.title.x = element_text(size = rel(1.25)),
                  axis.title.y = element_text(size = rel(1.25)),
                  plot.caption = element_text(face = "italic", size = 9),
                  legend.position = "bottom",
                  legend.background = element_rect(fill = "transparent"), 
                  panel.background = element_rect(fill = "transparent"),
                  panel.border = element_blank(),
                  axis.line = element_line(color="black"),
                  panel.grid.major.x = element_blank(),
                  panel.grid.minor.x = element_blank(),
                  panel.grid.minor.y = element_blank(),
                  panel.grid.major.y = element_line(color = "grey90", linetype = 2)
            ))
<<<<<<< HEAD

Last updated: 2021-02-19 14:40:16

=======

Last updated: 2021-02-18 15:14:59

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

Disclaimer

The Tennessee state information is simply a visual summary of the data directly downloaded from TN state Dept of Health, by date reported (rather than date of test). I do not attempt to model, make predictions, or give advice on how to use this data. If you want that sort of thing, go to the VUMC Health Policy site, which has developed some very good local models and Advisory Memos for Tennessee.

For the most up-to-date numbers visit the TN state Dept of Health or local health department sites. In particular, the Metro. Nashville site has a very nice summary- and note that their report of daily cases usually precedes the State Health Dept, due to the verification process and assignment of county of residence at the state.

Other useful sites for following trends in Tennessee and nationally:

TN Dept of Health Data

These data are taken directly from the TN state Dept of Health. I am using these data, by date of report, rather than date of test, for consistency.

# READ IN AND format data ####
tndoh_daily <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-Daily-Case-Info.xlsx"))
names(tndoh_daily) <- tolower(names(tndoh_daily))

# remove the new_test results which our bulk reported or adjusted
tndoh_daily <-
  tndoh_daily |> 
  mutate(date = lubridate::ymd(date)) |> 
  arrange(date) 
  
tndoh_daily[tndoh_daily$date %in% lubridate::ymd(c("2020-07-13", "2020-07-11", "2020-06-12")), "new_tests"] <- NA_real_

# STATE POPULATION DATA ####
us_data <-  readRDS(here::here("data_rds/us_data.rds"))
# state_pop <- readxl::read_excel(here::here("state_populations.xlsx"))
# state_pop$region <- gsub("\\.", "", state_pop$region)
# state_pop <-
#   state_pop |>
#   pivot_longer(cols = -c(region, census_2010),
#                names_to = "pop_year",
#                names_prefix = "pop_",
#                values_to = "pop_estimate")
# saveRDS(state_pop, here::here("data_rds/state_pop.rds"))

state_pop <- readRDS( here::here("data_rds/state_pop.rds"))
state_pop_2019 <-
  state_pop |> 
  filter(pop_year == 2019) |> 
  filter(region %in% state.name)
TN_pop2019 <- state_pop_2019$pop_estimate[state_pop_2019$region == "Tennessee"]

# vaccine state info
tndoh_vax_state <- readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_STATE_SUMMARY.XLSX"))
names(tndoh_vax_state) <- tolower(names(tndoh_vax_state))
tndoh_vax_state <-
  tndoh_vax_state %>% 
  mutate(date = lubridate::ymd(date),
         vaccine_count_pct = vaccine_count/TN_pop2019*100,
         recipient_count_pct = recipient_count/TN_pop2019*100,
         recip_fully_vacc_pct = recip_fully_vacc/TN_pop2019*100) 

# create annotations
tndoh_daily <-
  left_join(tndoh_daily, tndoh_vax_state,
            by = "date") |> 
  mutate(pos_tests_new = pos_tests -lag(pos_tests, 1), 
         pos_pct =   pos_tests/total_tests,
         pos_pct_today = pos_tests_new/new_tests,
         hospitalized_tot = replace_na(total_hosp, 0),
         hospitalized_new = total_hosp - lag(total_hosp),
         new_cases_100k = new_cases/TN_pop2019*100000,
         new_cases_7davg = RcppRoll::roll_meanr(new_cases, 7),
         new_cases_7davg_100k = new_cases_7davg/TN_pop2019*100000,
         new_deaths_7davg = RcppRoll::roll_meanr(new_deaths, 7),
         new_hosp_7davg = RcppRoll::roll_meanr(new_hosp, 7),
         new_tests_7davg = RcppRoll::roll_meanr(new_tests, 7),
         pos_pct_7davg = RcppRoll::roll_meanr(pos_tests_new, 7)/RcppRoll::roll_meanr(new_tests, 7),
         tt_daily = sprintf("<tr><td><b>%s</b></td></tr>
                     <tr><td><b>New Cases: %s</b></td></tr>
                     <tr><td><b>New Cases/100k: %s</b></td></tr>
                     <tr><td><b>New Hospitalizations: %s</b></td></tr>
                     <tr><td><b>New Deaths: %s</b></td></tr>
                     <tr><td><b>Test+ %% Today: %1.1f%%</b></td></tr>
                     <tr><td>New Cases 7d Avg/100k: %s</td></tr>
                     <tr><td>Deaths 7d Avg: %1.1f</td></tr>
                     <tr><td>Test+ %% 7d Avg: %1.1f%%</td></tr>
                     <tr><td>Case Fatality Ratio: %1.2f%%</td></tr>
                     <tr><td>New Vaccinations: %s</td></tr>
                     <tr><td>Rec 1st Vaccination %%Pop: %1.2f%%</td></tr>",
                     paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date), ", ", lubridate::year(date), " (", lubridate::wday(date, label = T), ")"), 
                     comma(new_cases, digits = 0), 
                     comma(new_cases_100k, digits = 0), 
                     comma(new_hosp, 0), 
                     comma(new_deaths, digits = 0), 
                     pos_pct_today*100, 
                     comma(new_cases_7davg_100k, 1),
                     comma(new_deaths_7davg, digits = 0), 
                     pos_pct_7davg*100, 
                     total_deaths/total_cases*100,
                     comma(new_vaccine_count, 0),
                     recipient_count_pct) ) 


tndoh_daily_limited <-
  tndoh_daily |>
  select(date, total_cases, total_hosp, total_deaths) |> 
  pivot_longer(cols = c(total_cases, total_hosp, total_deaths),
               names_to = "measure",
               values_to = "count") |> 
  mutate(var.f = factor(measure, levels = c("total_cases", "total_hosp", "total_deaths"),
                        labels = c("Cases", "Hospitalizations", "Deaths")))  
title_tnsdc <- cowplot::ggdraw() + 
  cowplot::draw_label(paste0("COVID-19 in Tennessee (as of ", max(tndoh_daily$date), ")") , fontface='bold')

tndoh_daily_long <-
  tndoh_daily |> 
  select(date, new_cases, new_deaths, new_hosp) |>  
  pivot_longer(cols = c(new_cases, new_deaths, new_hosp),
               names_to = "var", 
               values_to = "value") |> 
  mutate(var.f = factor(var, levels = c("new_cases", "new_deaths", "new_hosp"), 
                        labels = c("Cases", "Hospitalizations", "Deaths")))

label_pos_pct <-
    with(tndoh_daily,
         sprintf("Overall Positive: %1.1f%%\nToday's Positive: %1.1f%%", 
                 pos_pct[which(date == max(date))]*100,
                 pos_pct_today[which(date == max(date))]*100))



# AGE DATA ####
tndoh_age <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-Age.xlsx"))
names(tndoh_age) <-  tolower(names(tndoh_age))
tndoh_age <-
  tndoh_age |>
  filter(age_range != "Pending") |> 
  mutate(date = lubridate::ymd(date),
         age_range = gsub(" years", "", age_range),
         age.f = factor(age_range),
         deaths_case_rate = case_when(age.f == "Pending" ~ NA_real_,
                                      TRUE ~ ar_totaldeaths/ar_casecount*100))

# TN age groups 2019 estimates(taken from https://www.census.gov/data/tables/time-series/demo/popest/2010s-state-detail.html)
tn_age_groups <- 
  tibble::tibble(age.f = factor(c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"),
                                    levels = c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")),
                 populationtn_byage_2019 = c(821984L, 855582L, 942211L, 887898L, 850125L, 893578L, 808164L, 516865L, 252767L)) 

tndoh_age <- 
  left_join(tndoh_age,
          tn_age_groups,
          by = "age.f") |> 
  group_by(age.f) |> 
  mutate(ar_new_cases_100k = new_arcases/populationtn_byage_2019*100000,
         ar_new_cases_100k_7davg = round(RcppRoll::roll_meanr(new_arcases, 7)/populationtn_byage_2019*100000, 1),
         ar_totaldeaths_100k    = ar_totaldeaths/populationtn_byage_2019*100000,
         ar_new_deaths_100k = ar_newdeaths/populationtn_byage_2019*100000,
         ar_new_deaths_100k_7davg = round(RcppRoll::roll_meanr(ar_newdeaths, 7)/populationtn_byage_2019*100000, 1)) 



# VACCINE by age
 tndoh_vax_age <- readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_AGE_GROUPS.XLSX"),
                                      col_types = c("text", "date", rep("numeric", 6)))
names(tndoh_vax_age) <- tolower(names(tndoh_vax_age))
tndoh_vax_age <-
  tndoh_vax_age %>%
  filter(!grepl("pending", age_group, ignore.case = T )) |> 
  mutate(date = lubridate::ymd(date),
         age.f = factor(age_group,
                        levels = c("0-10", "11-20", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"), # 16-20 needs to match other age category
                        labels = c("0-10", "11-20", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")))
tndoh_age <-
  left_join(tndoh_age, tndoh_vax_age,
            by = c("date", "age.f")) %>%
  rowwise() %>%
  mutate(ar_totalvax_pct = vaccine_count/populationtn_byage_2019*100,
         ar_recipientvax_pct = recipient_count/populationtn_byage_2019*100) %>%
  ungroup() 

# most recent vaccination by age
vax_by_age_mostrecent <-
  tndoh_age %>% 
  filter(!is.na(vaccine_count)) %>% 
  filter(date == max(date, na.rm = T)) %>% 
    select(date, recipient_count, age.f, ar_recipientvax_pct) 
mostrecent_vax_agedate <- max(vax_by_age_mostrecent$date)  

# Deaths by age table ####
age_death_table <-
  tndoh_age %>%
  filter(date == max(date, na.rm = T)) %>% 
  select(age.f,populationtn_byage_2019, ar_casecount, ar_totaldeaths, ar_totaldeaths_100k) %>% 
  mutate(#ar_casecount = comma(ar_casecount, 0),
    ar_totaldeaths_pct = sprintf("%1.2f%%", ar_totaldeaths_100k/1000),
    ar_totaldeaths_100k = sprintf("%1.2f", ar_totaldeaths_100k)) %>% 
  left_join(vax_by_age_mostrecent, by = "age.f")


# COUNTY DATA ####
# tndoh_countynew <- readxl::read_excel(here::here("tn_data_center/TN_COVID19_CountyDaily.xlsx"))
TN_counties <- readRDS(here::here("data_rds/TN_counties.rds")) # see below for source

tndoh_county_coln <-   length(names(readxl::read_excel(here::here("data_tndoh/Public-Dataset-County-New.xlsx"), n_max = 2)))

tndoh_countynew <-  readxl::read_excel(here::here("data_tndoh/Public-Dataset-County-New.xlsx"),
                                       col_types = c("date", "text", rep("numeric", tndoh_county_coln-2))) 

names(tndoh_countynew) <-  tolower(names(tndoh_countynew))
tndoh_countynew <-
  tndoh_countynew |> 
  mutate(date = lubridate::ymd(date),
         pct_pos = total_cases/(total_tests)*100,
         county_lower = tolower(county)) |> 
  group_by(county) |>
  arrange(date) |> 
  mutate(new_cases_county_4davg = round(RcppRoll::roll_meanr(new_cases, 4), 1),
         new_cases_county_7davg = round(RcppRoll::roll_meanr(new_cases, 7), 1),
         new_deaths_county_7davg = round(RcppRoll::roll_meanr(new_deaths, 7), 1),
         new_testpos_pct   = round(new_pos_tests/new_tests*100, 1),
         new_testpos_pct_7davg   = round(RcppRoll::roll_meanr(new_pos_tests, 7)/RcppRoll::roll_meanr(new_tests, 7)*100, 1),
         new_hospitalized_county_7davg = round(RcppRoll::roll_meanr(new_hospitalized, 7), 1),
         new_cases_county_7davg_1wkdelta = (new_cases_county_7davg - lag(new_cases_county_7davg, 7))/lag(new_cases_county_7davg, 7) 
  ) |> 
  ungroup()


  
# TN_data <- tigris::counties(state = "TN", cb = TRUE)
# names(TN_data) <- tolower(names(TN_data))
# TN_counties <- sf::st_as_sf(TN_data) |> 
#     mutate(county_lower = tolower(as.character(name)))

# TN COUNTY POPULATION FROM US CENSUS 2019 ESTIMATES
# tn_county_pop <- read_csv(here::here("data_csv/TN 2019 county pop estimates.csv"))
# tn_county_pop <-
#   tn_county_pop |> 
#   mutate(county_lower = tolower(county)) |> 
#   select(-county)
# TN_counties <-
#   left_join(TN_counties,
#             tn_county_pop,
#             by = "county_lower")


# saveRDS(TN_counties, here::here("data_rds/TN_counties.rds"))
# TN_counties <- readRDS(here::here("data_rds/TN_counties.rds"))

tndoh_countynew <-
  left_join(tndoh_countynew,
            TN_counties,
            by = "county_lower") |> 
  mutate(total_cases_per100kpop = total_cases/population*100000,
         new_cases_per100kpop_7davg   = new_cases_county_7davg /population*100000,
         total_tests_pctpop = total_tests/population*100,
         new_tests_per100kpop   = new_tests/population*100000,
         total_deaths_per100kpop = total_deaths/population*100000) 

tndoh_countynew_calc <-
  tndoh_countynew |> 
    group_by(date) |> 
    filter(tolower(county) %in% TN_counties$county_lower) |> # exclude out of state and pending
    mutate(#new_testpos_pct = new_pos_tests/new_tests*100,
           new_cases_county_7davg_rank = rank(-new_cases_county_7davg),
           new_cases_per100kpop_7davg_rank = rank(-new_cases_per100kpop_7davg),
           new_cases_rank = rank(-new_cases)) |> 
    ungroup() |> 
    select(county, date, new_cases_county_7davg_rank, 
           new_cases_per100kpop_7davg_rank, new_cases_rank)
tndoh_countynew <-   
  left_join(tndoh_countynew,
            tndoh_countynew_calc,
            by = c("date", "county"))  

 

# DEMOGRAPHICS
tndoh_demographics <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-RaceEthSex.xlsx"))
names(tndoh_demographics) <-  tolower(names(tndoh_demographics))
tndoh_demographics$cat_detail <- gsub("/ ", "/", tndoh_demographics$cat_detail)
tndoh_demographics <-
  tndoh_demographics |> 
  filter(category == "RACE",
         date == max(date)) |> 
  rename(race = cat_detail,
         cases_byrace = cat_casecount) |> 
  mutate(date = lubridate::ymd(date),
         race = factor(race,
                       levels = c("White", "Black or African American", "Asian",  "Other/Multiracial", "Pending",  "American Indian or Alaska Native", "Native Hawaiian or Other Pacific Islander"),
                       labels = c("White", "Black", "Asian",  "Other", "Pending",  "Am.Indian/Alaskan", "Hawaiian/Pacific")))


# US DATA FROM THE COVIDTRACKING PROJECT ####
covus_currentdata <- readr::read_csv(here::here("data_covidtracking/daily.csv"))
covus_currentdata <-
  covus_currentdata |> 
  filter(state %in% state.abb,
         !state %in% c("AK", "HI") ) |> 
  mutate(date = lubridate::ymd(date))

# covus_currentdata <- readRDS(here::here("data_covidtracking/covus_current.rds"))

# us_data <- tigris::states(cb = TRUE)
# us_data <-   sf::st_as_sf(us_data)
# names(us_data) <- tolower(names(us_data))
# us_data <-
#   us_data |> 
#   filter(stusps %in% state.abb)
# saveRDS(us_data, here::here("data_rds/us_data.rds"))
# us_data <-  readRDS(here::here("data_rds/us_data.rds"))
# 
# 
# 
# # state_pop <- readxl::read_excel(here::here("state_populations.xlsx"))
# # state_pop$region <- gsub("\\.", "", state_pop$region)
# # state_pop <-
# #   state_pop |>
# #   pivot_longer(cols = -c(region, census_2010),
# #                names_to = "pop_year",
# #                names_prefix = "pop_",
# #                values_to = "pop_estimate")
# # saveRDS(state_pop, here::here("data_rds/state_pop.rds"))
# 
# state_pop <- readRDS( here::here("data_rds/state_pop.rds"))
# state_pop_2019 <-
#   state_pop |> 
#   filter(pop_year == 2019) |> 
#   filter(region %in% state.name)

covus_currentdata <-
  covus_currentdata |> 
  left_join(us_data,
            by = c("state" = "stusps")) |> 
    left_join(state_pop_2019, by = c("name" = "region")) |> 
  mutate(test_density_pct = total/pop_estimate*100,
         case_density_pct = positive/pop_estimate*100)  

covus_currentdata <-
  covus_currentdata |>
  group_by(state) |> 
  arrange(date) |> 
  mutate(new_cases_state_7davg = round(RcppRoll::roll_meanr(positiveIncrease, 7), 1),
         new_deaths_state_7davg = round(RcppRoll::roll_meanr(deathIncrease, 7), 1),
         new_deaths_state_7davg_100k = new_deaths_state_7davg/pop_estimate*100000,
         new_testpos_pct   = round(positiveIncrease/totalTestResultsIncrease*100, 1),
         new_testpos_pct_7davg   = round(RcppRoll::roll_meanr(positiveIncrease, 7)/RcppRoll::roll_meanr(totalTestResultsIncrease, 7)*100, 1),
         new_cases_state_7davg_100k = new_cases_state_7davg/pop_estimate*100000
  ) |> 
  ungroup()

covus_currentdata <-
  covus_currentdata |>
  group_by(date) |> 
  mutate(rank_newcases7davg = rank(-new_cases_state_7davg_100k)) |> 
  ungroup()



state_tbl <- 
  tibble(state_name = state.name, 
         state = state.abb)


covus_currentdata_2163 <-
  covus_currentdata |> 
  mutate(geometry = st_transform(geometry, 2163)) |> 
    left_join(state_tbl, by = "state") |> 
  mutate(tt = sprintf("<tr><td><b><i>%s (%s)</i></b></td></tr>
                     <tr><td>New Cases/100k (7d Avg): %s</td></tr>
                     <tr><td>Test %% Positive (7d Avg): %1.1f%%</td></tr>
                     <tr><td>State Rank (Cases/100k 7d Avg): %s</td></tr>
                     <tr><td>Total Cases: %s</td></tr>
                     <tr><td>Total Deaths: %s</td></tr>
                     <tr><td>New Deaths/d/100k (7d Avg): %1.2f</td></tr>
                     <tr><td>CFR: %1.2f%%</td></tr>
                     <tr><td>Case Density: %1.1f%% of pop.</td></tr>
                     <tr><td>Testing Density: %1.1f%% of pop.</td></tr>", 
                     toupper(state_name),  paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date)),
                     comma(new_cases_state_7davg_100k, digits = 1), new_testpos_pct_7davg,  rank_newcases7davg, comma(positive, digits = 0),
                     comma(death, digits = 0), new_deaths_state_7davg_100k, death/positive*100, case_density_pct, test_density_pct)) 



# county
tndoh_vax_county <- readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_COUNTY_SUMMARY.XLSX"),
                       col_types = c("date", "text", rep("numeric", 6)))
    

names(tndoh_vax_county) <- tolower(names(tndoh_vax_county))
tndoh_vax_county <-
  tndoh_vax_county %>% 
  mutate(date = lubridate::ymd(date),
         county_lower = tolower(county)) %>% 
    left_join(TN_counties,
              by = "county_lower") %>% 
    rowwise() %>%  
    mutate(vaccine_count_pct = vaccine_count/population*100,
           recipient_count_pct = recipient_count/population*100) %>% 
    ungroup() %>% 
    filter(!is.na(statefp))

Testing Update- Interactive

interactive version- mouse/hover over bars for details

pl_new_tnsdc <-
  tndoh_daily |>
  pivot_longer(cols = c(new_cases, new_hosp, new_deaths),
               names_to = "var", 
               values_to = "value") |> 
  ggplot(aes(x = date)) |> \(x) x +
  ggiraph::geom_col_interactive(aes(y = value, fill = var, tooltip = tt_daily, data_id = date), alpha = 0.5,
                               position = "identity") +
  geom_line(aes(y = new_cases_7davg),color = "steelblue", size = 1) +
  labs(x=NULL,
       y="Results by day reported",
       subtitle = "Results by Day, New (Total)",
       fill = NULL) +
  scale_fill_manual(values = c("steelblue",
                               "red",
                               "orange"),
                    labels = c(paste0("Cases: ", comma(tndoh_daily$new_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
                                      " (", comma(tndoh_daily$total_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
                               paste0("Deaths: ", comma(tndoh_daily$new_deaths[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
                                      " (", comma(tndoh_daily$total_deaths[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
                               paste0("Hospitalizations: ", comma(tndoh_daily$new_hosp[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
                                      " (", comma(tndoh_daily$total_hosp[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")")
                    )) +
  scale_y_continuous(expand = c(0,0),
                     sec.axis = sec_axis( trans=~./TN_pop2019*10^5, name="Cases/100k Population")) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b" ) +
  theme(
    legend.position = c(.01, 1), 
    legend.justification = c(0, 1),
    legend.text = element_text(size = 10))

hover_css <- "background-color:orange;stroke:black;fill:orange;"
tooltip_css <- " background-color:rgba(255, 255, 255, 0.9);color:black;stroke:black;opacity:0.1;border:2px solid black;border-radius: 8px;font-family:Arial;font-size: 0.875em;"

ggiraph::girafe(ggobj = pl_new_tnsdc,
                width_svg = 8,
                options = list(
                  offx = 40, offy = 40,
                  opts_hover(css = hover_css),
                  opts_tooltip(css = tooltip_css)
                ) )
<<<<<<< HEAD
=======
>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

Today’s Overview

# PLOTS ####
pl_tests_new <-
  tndoh_daily |>
  ggplot(aes(x=date)) |> \(x) x +
  geom_col(aes(y = new_tests, fill = "Tests"), alpha = 0.5) +
  geom_col(aes(y = new_cases, fill = "Cases"), alpha = 0.7) +
  labs(x=NULL,
       fill = NULL,
       y="Tests by day reported",
       subtitle = "Tests by Day, New (Total)") +
  scale_fill_manual(values = c("Cases" = "steelblue",
                               "Tests" = "grey"),
                      labels = c(paste0("Cases:", comma(tndoh_daily$new_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
                                        " (", comma(tndoh_daily$total_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
                                 paste0("Tests:", comma(tndoh_daily$new_tests[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), 
                                        " (", comma(tndoh_daily$total_tests[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"))) +
  scale_y_continuous(expand = c(0,0),
                     sec.axis = sec_axis( trans=~./TN_pop2019*10^5, name="per 100k Population")) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b" ) +
  annotate(geom = "text", x=min(tndoh_daily$date)+1, y = 0.65*max(tndoh_daily$new_tests, na.rm = T), hjust = 0, 
           label = label_pos_pct, size = 3.5) +
  theme(#axis.text.x =  element_text(angle = 45, hjust = 1),
        legend.position = c(.01, 1), 
        legend.justification = c(0, 1),
        legend.text = element_text(size = 10))


# age ####

pl_age <-
  tndoh_age |>
  filter(date == max(date)) |> 
  ggplot() |> \(x) x +
# pl_age <-
#   pl_age +
  geom_blank(data = tndoh_age[tndoh_age$ar_casecount == max(tndoh_age$ar_casecount), ], 
             aes(x=age.f, y=ar_casecount*1.1)) +
  geom_bar(aes(age.f, ar_casecount), alpha = 0.5, fill = "steelblue", stat = "identity", ) +
  geom_bar(aes(age.f, ar_totaldeaths), alpha = 0.5, fill = "red", stat = "identity", ) +
  geom_text(aes(age.f, ar_casecount*1.02, vjust =0, label = comma(ar_casecount, digits = 0)),
            color = "steelblue", size = 3) +
    geom_text(aes(age.f, ar_totaldeaths*1.01, vjust =0,  label = comma(ar_totaldeaths, digits = 0)),
              color = "red", size = 3) +
    labs(x=NULL,
       y="Total Case Count",
       subtitle = "Total Cases by Age") +
  scale_y_continuous(expand = c(0,0)) +
  theme(axis.text.x =  element_text(angle = 45, hjust = 1)) 

# multi-axis plot
pl_age_death <-
  tndoh_age |> 
    filter(date == max(date)) |> 
  ggplot(aes(age.f, deaths_case_rate, group = date)) |> \(x) x +
# pl_age_death <-
#   pl_age_death +
  geom_point(color = "red") +
  geom_line(color = "red")

# Dual axis 
ggp <- ggplot_build(pl_age)
ggp2 <- ggplot_build(pl_age_death)

ylim.prim <- ggp$layout$panel_params[[1]]$y.range   # in this example, precipitation
ylim.sec <- ggp2$layout$panel_params[[1]]$y.range    # in this example, temperatureb <- diff(ylim.prim)/diff(ylim.sec)
b <- diff(ylim.prim)/diff(ylim.sec)
a <- b*(ylim.prim[1] - ylim.sec[1])

pl_age_death_final <-
  tndoh_age |>
  filter(date == max(date)) |> 
  ggplot(aes(x=age.f)) |> \(x) x +
# pl_age_death_final <-
#   pl_age_death_final +
  geom_blank(data = tndoh_age[tndoh_age$ar_casecount == max(tndoh_age$ar_casecount), ], 
             aes(x=age.f, y=ar_casecount*1.1)) +
  geom_bar(aes(age.f, ar_casecount), alpha = 0.5, fill = "steelblue", stat = "identity", ) +
  geom_bar(aes(age.f, ar_totaldeaths), fill = "red", stat = "identity", ) +
  geom_line(aes(y = a + deaths_case_rate*b, group = date), color = "red", size = 1.5) +
  geom_point(aes(y = a + deaths_case_rate*b, group = date), color = "red", size = 2, shape = 21, fill = "white") +
  geom_text(aes(age.f, ar_casecount*1.02, vjust =0, label = comma(ar_casecount, digits = 0)),
            color = "steelblue", size = 3) +
    geom_text(aes(age.f, ar_totaldeaths*1.01, vjust =0,  label = comma(ar_totaldeaths, digits = 0)),
              color = "red", size = 3) +
    labs(x=NULL,
       y="Total Case Count",
       subtitle = "Total Cases by Age") +
  scale_y_continuous(expand = c(0,0),
                     sec.axis = sec_axis(~ (. - a)/b, name = "Crude Case Fatality (%)")) +
  theme(axis.text.x =  element_text(angle = 45, hjust = 1),
        axis.title.y.right = element_text(color = "red"),
        axis.text.y.right = element_text(color = "red")) 


race_total <- sum(tndoh_demographics$cases_byrace[tndoh_demographics$race != "Pending"])

# AGE HEATMAP
max_newcaserate_byage <- max(tndoh_age$ar_new_cases_100k_7davg, na.rm = T)
ar_newcase_heatmap_100k <-
  tndoh_age |> 
  filter(date >= lubridate::ymd("2020-04-01")) |> 
  ggplot(aes(x = date,
             y = age.f)) |> \(x) x +
# ar_newcase_heatmap_100k <-
#   ar_newcase_heatmap_100k +
  geom_tile(aes(fill = ar_new_cases_100k_7davg,
                color = ar_new_cases_100k_7davg)) +
  viridis::scale_fill_viridis(labels = scales::comma,
                              limits = c(0, max_newcaserate_byage), oob = scales::squish,
                              option = "plasma", na.value="black") +
  viridis::scale_color_viridis(labels = scales::comma,
                               limits = c(0, max_newcaserate_byage), oob = scales::squish,
                               option = "plasma", na.value="black") +
  labs(fill = "New Cases/100k/d\n(7dAvg)",
       color  = "New Cases/100k/d\n(7dAvg)",
       x = NULL, y = NULL,
       subtitle = "New Case Rate, by Age (population adjusted)") +
  scale_x_date(date_breaks = "1 month", date_labels = "%b", expand = c(0, 0)) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 12),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.background = element_blank(), 
        panel.border = element_blank(),
        legend.title = element_text(size = 10),
        legend.text =  element_text(size = 8))

# fig 1 ####  
pl_TOP <- cowplot::plot_grid(pl_new_tnsdc, pl_tests_new, 
                             pl_age_death_final, ar_newcase_heatmap_100k,
                             align = "h", axis = "bt",
                             nrow = 2)
cowplot::plot_grid(title_tnsdc, pl_TOP, nrow = 2, rel_heights = c(.1, 1))
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 
# ggsave(here::here("figures/pl_tnmap.jpg"), dpi = 600, width = 12, height = 4, units = "in")

SUMMARY BY AGE In Tennessee, 1.68% of residents >80 and 0.65% of those age 71-80 have died of COVID-19, to date.

Age Group Population (n) Cases Cases
(% of Total) 		Deaths Deaths
(% of Total) 		Deaths/100k Deaths (% in Age) Given 1st Vaccine (% in Age)
0-10 821984 <<<<<<< HEAD 40669 ======= 40596 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 5% 4 0% 0.49 0.00% NA%
11-20 855582 <<<<<<< HEAD 96383 ======= 96168 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 13% 4 0% 0.47 0.00% 0.5%
21-30 942211 <<<<<<< HEAD 138036 ======= 137842 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 18% 45 0% 4.78 0.00% 5.6%
31-40 887898 <<<<<<< HEAD 118184 ======= 117988 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 16% 111 1% 12.50 0.01% 7.6%
41-50 850125 <<<<<<< HEAD 113616 ======= 113417 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 15% 346 3% 40.70 0.04% 8.6%
51-60 893578 <<<<<<< HEAD 107155 ======= 106943 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 14% 915 8% 102.40 0.10% <<<<<<< HEAD 9.1% ======= 9.0% >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e
61-70 808164 <<<<<<< HEAD 77470 ======= 77314 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 10% <<<<<<< HEAD 2013 ======= 2011 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 18% <<<<<<< HEAD 249.08 ======= 248.84 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 0.25% <<<<<<< HEAD 12.3% ======= 12.0% >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e
71-80 516865 <<<<<<< HEAD 46031 ======= 45954 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 6% <<<<<<< HEAD 3382 ======= 3380 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 31% <<<<<<< HEAD 654.33 ======= 653.94 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 0.65% <<<<<<< HEAD 44.9% ======= 44.3% >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e
81+ 252767 <<<<<<< HEAD 24078 ======= 24026 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 3% <<<<<<< HEAD 4244 ======= 4241 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 38% <<<<<<< HEAD 1679.02 ======= 1677.83 >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 1.68% <<<<<<< HEAD 48.2% ======= 47.9% >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e
<<<<<<< HEAD Population within each age group taken from US Census estimates for TN, 2019
Vaccination data last updated on 2021-02-19. ======= Population within each age group taken from US Census estimates for TN, 2019
Vaccination data last updated on 2021-02-18. >>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

NOTES:

  • “new_test” numbers reported on “2020-07-13”, “2020-07-11”, and “2020-06-12” are not representative of the time trend (reclassification of how tests were reported, dumped on a a few dates) and are omitted
  • the reporting format changed to include “probable” cases on 6/12/2020. A total of 202 cases were classified as “probable” up to that date, and they were all included in the new case count on 6/12/2020.
  • Test reporting Note: Beginning on 6/12/2020 the number of all tests reported includes all tests, including repeat tests in the same person; an unclear number of tests were added on 6/12/2020 to account for this change in reporting.
  • COVID-19 data reported for June 29, 2020 reflect results from 2 days- due to an unplanned shutdown of the state surveillance system during the weekend.
  • Age strata populations are taken from the US Census estimate for TN, year 2019 estimate (the most recent available)

Vaccinations

A bit of good news. Vaccination data is not updated every day- most recent data provided.

TN now provides detail on number of recipients and number of vaccinations.

latest_vax_date <- max(tndoh_vax_county$date, na.rm = T)
tndoh_vax_county %>% 
  filter(date == latest_vax_date,
         !is.na(statefp)) %>% 
    ggplot() |> \(x) x +
    geom_sf(data = TN_counties) +
    geom_sf(aes(geometry = geometry,
                fill = recipient_count_pct),
            color = "grey30",
            inherit.aes = F) +
    geom_sf_text(aes(label=sprintf("%0.1f%%", recipient_count_pct), 
                     # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos), 
                     geometry = geometry, 
                     color = cut(recipient_count_pct, 2)), 
                 size = 2.5) + 
    viridis::scale_fill_viridis(option = "C", direction = 1) +
    scale_color_manual(values = c("white", "black"), guide = F) +
    labs(fill = "%Pop.:",
         title = paste0("Vaccination Coverage (latest update on ", latest_vax_date, ")"),
         subtitle = sprintf("First Vaccinations: %s (%1.1f%% of TN Population)\nCompleted Vacc. : %s (%1.1f%% of Pop.)\nTotal Vaccinations: %s",
                            comma(max(tndoh_vax_state$recipient_count, na.rm = T), 0),
                            max(tndoh_vax_state$recipient_count_pct, na.rm = T),
                            comma(max(tndoh_vax_state$recip_fully_vacc, na.rm = T), 0),
                            max(tndoh_vax_state$recip_fully_vacc_pct, na.rm = T),
                            comma(max(tndoh_vax_state$vaccine_count, na.rm = T), 0))
         ) +
    theme_void() +
    theme(legend.position = "left")
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

Rolling 7-Day Averages & Top 10 Most Cases

note: The “New Hospitalizations” in this chart, as reported by the state, may not accurately reflect actual hospitalizations. Right y-axis is population adjusted (per 100k for TN). Ignore for “Test Positivity.”

Counties in the top 10 highest new case rate over the past 7 days; 7-day average trend-line for daily new cases

variable.labs <- c( "Cases", "Hospitalizations", "Deaths", "Tests", "Test Positivity (%)", "Vaccinations")
names(variable.labs) <- c( "new_cases", "new_hosp", "new_deaths", "new_tests", "pos_pct_today", "new_vaccine_count")

plot_state_trends <- 
  tndoh_daily |> 
  mutate(pos_pct_today = pos_pct_today*100) |> 
  # select(date, new_cases_7davg)
  # select(date, new_cases, new_deaths, new_tests, pos_pct_today) |> 
  select(date, new_cases, new_hosp, new_deaths, new_tests, pos_pct_today, new_vaccine_count) |>
  pivot_longer(cols = -date,
               names_to = "variable",
               values_to = "value") |> 
  group_by(variable) |> 
  arrange(date) |> 
  mutate(value_7davg = RcppRoll::roll_meanr(value, 7)) |> 
  ggplot() |> \(x) x +
# plot_state_trends <- 
#   plot_state_trends +
  geom_line(aes(date, value), size = 0.6, color = "grey") +
  geom_line(aes(date, value_7davg), size = 1, color = "red") +
  scale_y_continuous(sec.axis = sec_axis( trans=~./TN_pop2019*10^5, name="per 100k Population")) +
  labs(x=NULL,
       y="Rolling 7-day Average",
       title = "TN State Trends") +
  facet_wrap(~variable, scales = "free_y", nrow = 2,
             labeller = labeller(variable = variable.labs)) +
  theme(strip.text = element_text(size = 14))


# MIDTN counties plot ####
counties_top10 <-
  tndoh_countynew |>  #select(new_cases_county_7davg)
  ungroup() |> 
  filter(date == max(date, na.rm = T),
         name %in% TN_counties$name) |> # exclude the pending and out-of-state
  top_n(10, wt = new_cases_county_7davg) |> 
  pull(name)

plot_top10_counties <- 
  tndoh_countynew |> 
  filter(name %in% counties_top10) |>
  ggplot() |> \(x) x +
# plot_top10_counties <- 
#   plot_top10_counties +
  # geom_bar(stat = "identity") +
  geom_col(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
  geom_col(aes(date, new_deaths), fill = "red") +
  geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
  labs(x=NULL,
       y="New Results by day",
       fill = NULL,
       alpha = NULL,
       title = "Top 10 New Cases by Day (7d avg.)") +
  scale_y_continuous(expand = c(0,0)) +
  ylim(0,NA) +
  scale_x_date(date_breaks = "3 month", date_labels = "%b" ) +
  theme(#axis.text.x =  element_text(angle = 90, hjust = 1),
        legend.position = "none") +
  facet_wrap(~name, nrow = 2, 
             scales = "free_y")+
  theme(strip.text = element_text(size = 14))

cowplot::plot_grid(plot_state_trends,
                   plot_top10_counties,
                   nrow = 2, rel_heights = c(1.5,2))
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

Top 10: Largest Increase in Cases

Counties with the top 10 greatest weekly increase in New Cases, by % Change of 7-day average versus 7 days prior (limited to counties with average of >5 daily cases).

# MIDTN counties plot ####
counties_top10_increase <-
  tndoh_countynew |>  
  ungroup() |> 
  filter(date == max(date, na.rm = T),
       new_cases_county_7davg_1wkdelta > 0,
         new_cases_county_7davg >= 5,
         name %in% TN_counties$name) |> # exclude the pending and out-of-state
  top_n(n = 10, wt = new_cases_county_7davg_1wkdelta) |> 
  pull(name)

plot_top10_increase <-
  tndoh_countynew |> 
  filter(name %in% counties_top10_increase) |>
  ggplot() |> \(x) x +
# plot_top10_increase + 
  # geom_bar(stat = "identity") +
  geom_col(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
  geom_col(aes(date, new_deaths), fill = "red") +
  geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
  labs(x=NULL,
       y="New Results by day",
       fill = NULL,
       alpha = NULL,
       subtitle = "New Cases by Day: 10 Counties with the largest 1-week increase") +
  scale_y_continuous(expand = c(0,0)) +
  ylim(0,NA) +
  scale_x_date(date_breaks = "3 month", date_labels = "%b" ) +
  theme(#axis.text.x =  element_text(angle = 90, hjust = 1),
        legend.position = "none") +
  facet_wrap(~name, nrow = 2, 
             scales = "free_y")+
  theme(strip.text = element_text(size = 14))

plot_top10_increase
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

Best 10: Largest Decrease in Cases

Counties with the top 10 greatest weekly Decrease in New Cases, by % Change of 7-day average versus 7 days prior (limited to counties with average of >5 daily cases).

Only counties with a decline in case averages are shown (may be <10)

# MIDTN counties plot ####
counties_top10_decrease <-
  tndoh_countynew |>  
  ungroup() |> 
  filter(date == max(date, na.rm = T),
       new_cases_county_7davg_1wkdelta < 0,
           new_cases_county_7davg >= 5,
         name %in% TN_counties$name) |> # exclude the pending and out-of-state
  top_n(n = -10, wt = new_cases_county_7davg_1wkdelta) |> 
  pull(name)

pl_best10 <-
 tndoh_countynew |> 
  filter(name %in% counties_top10_decrease) |>
  ggplot() |> \(x) x +
# pl_best10 <-
#   pl_best10 +
  # geom_bar(stat = "identity") +
  geom_col(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
  geom_col(aes(date, new_deaths), fill = "red") +
  geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
  labs(x=NULL,
       y="New Results by day",
       fill = NULL,
       alpha = NULL,
       subtitle = "New Cases by Day: 10 Counties with the largest 1-week Decrease") +
  scale_y_continuous(expand = c(0,0)) +
  ylim(0,NA) +
  scale_x_date(date_breaks = "3 month", date_labels = "%b" ) +
  theme(#axis.text.x =  element_text(angle = 90, hjust = 1),
        legend.position = "none",
        strip.text = element_text(size = 14)) +
  facet_wrap(~name, nrow = 2, 
             scales = "free_y")
pl_best10
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 
# cat("There were no counties with a decreased 7d average")


# ifelse(length(counties_top10_decrease) >= 1,
#        # IF >=1 with negative trend
#        plot(pl_best10),
#        # IF No counties with negative trend
#        "There are no counties with a decreasing case trend")

Interactive Tennessee Map

County color indicates the New Cases per day per county population (7 day average)

# make county data to list
county_plots <-
  tndoh_countynew |>
  filter(!is.na(date),
         date >= lubridate::ymd("2020-03-15")) |> 
  group_by(name) |> 
  nest() |> 
  arrange(name) |> 
  mutate(pl_new = map2(data, name,
                       ~ggplot(data = .x) +
                         geom_col(aes(date, new_cases, fill = "Cases"), alpha = 0.5) +
                         geom_col(aes(date, new_deaths, fill = "Deaths")) +
                         geom_line(aes(date, new_cases_county_7davg, fill = "Case 7-Day Avg"), color = "steelblue", size = 1) +
                         # geom_col(aes(date, new_hospitalized, fill = "Hospitalized")) +
                         labs(x=NULL,
                              y="New Results by day",
                              title = glue::glue("{.y} County"),
                              fill = NULL,
                              alpha = NULL,
                              subtitle = "New Results by Date Reported") + 
                         scale_fill_manual(values = c("Cases" = "steelblue",
                                                      "Case 7-Day Avg" = "steelblue",
                                                      "Deaths" = "red"),
                                           # "Hospitalized" = "orange"),
                                           labels = c(paste0("Total Cases: ", comma(.x$total_cases[.x$date == max(.x$date)], digits = 0)),
                                                      paste0("7d Avg: ", comma(.x$new_cases_county_7davg[.x$date == max(.x$date)], digits = 0)),
                                                      paste0("Total Deaths: ", comma(.x$total_deaths[.x$date == max(.x$date)], digits = 0))
                                                      # paste0("Hospitalized: ", .x$total_hospitalized[.x$date == max(.x$date)])
                                           )
                         ) +
                         scale_y_continuous(expand = c(0,0)) +
                         ylim(0,NA) +
                         scale_x_date(date_breaks = "1 month", date_labels = "%m/%d" ) +
                         theme(axis.text.x =  element_text(angle = 45, hjust = 1),
                               legend.position = "bottom")
  ))

county_data_byday_now <-
  left_join(TN_counties, 
            {tndoh_countynew |> 
                filter(date == max(date, na.rm = T)) |>  
                select(name, total_cases, new_cases_county_7davg, new_cases_per100kpop_7davg, new_cases_per100kpop_7davg_rank,
                       new_testpos_pct, new_testpos_pct_7davg, new_tests, new_pos_tests,
                       new_cases_rank, new_cases, new_cases_rank,
                       total_deaths, total_hospitalized)}) |> 
  mutate(cases = replace_na(total_cases, 0),
         deaths = replace_na(total_deaths, 0)) |> 
  left_join({county_plots |> select(name, pl_new)},
            by = "name")


lbls <-
  sprintf("<strong><i>%s Co.</i></strong><br/>
          Total Cases: %s<br/>
          New Cases/d/100k, 7dAvg (rank): %1.1f (%s/%s)<br/>
          New Cases Today (rank): %s (%s/%s)<br/>
          Test +%% Today: %1.1f%% (%s/%s)<br/>
          Test +%% (7dAvg): %1.1f%%<br/>
          Total Deaths: %s<br/>
          <i>Click to view time-course</i>", 
          toupper(county_data_byday_now$name), 
          comma(county_data_byday_now$total_cases, digits = 0), 
          county_data_byday_now$new_cases_per100kpop_7davg, county_data_byday_now$new_cases_per100kpop_7davg_rank, max(county_data_byday_now$new_cases_per100kpop_7davg_rank), 
          comma(county_data_byday_now$new_cases, digits = 0),county_data_byday_now$new_cases_rank, max(county_data_byday_now$new_cases_rank), 
          county_data_byday_now$new_testpos_pct, county_data_byday_now$new_pos_tests, county_data_byday_now$new_tests, 
          county_data_byday_now$new_testpos_pct_7davg,
          county_data_byday_now$total_deaths) |> 
  lapply(htmltools::HTML)


# pal <- colorQuantile("viridis", domain = county_data_byday_now$new_cases_per100kpop_7davg, n = 10, reverse = FALSE)
pal <- colorNumeric("viridis", domain = county_data_byday_now$new_cases_per100kpop_7davg, reverse = FALSE)


leaflet() |> 
  addProviderTiles(providers$Stamen.TonerLite) |>
  addPolygons(data = TN_counties, # borders of all counties
              color = "grey",
              fill = NA,
              weight = 1) |>
  addPolygons(data = county_data_byday_now,
              group = "name",
              fillColor = ~pal(county_data_byday_now$new_cases_per100kpop_7davg),
              fillOpacity = 0.8,
              weight = 1,
              highlight = highlightOptions(
                weight = 5,
                color = "red",
                bringToFront = TRUE),
              label = lbls,
              labelOptions = labelOptions(
                style = list("font-weight" = "normal", padding = "3px 8px"),
                textsize = "15px",
                direction = "auto")) |>
  addLegend(pal = pal, 
            values = county_data_byday_now$new_cases_per100kpop_7davg,
            title = "New Cases/100k (7d Avg)",
            opacity = 0.9,
            position = "bottomright") |> 
  leafpop::addPopupGraphs(county_data_byday_now$pl_new,
                          group = "name",
                          width = 400, height = 200)
<<<<<<< HEAD
=======
>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

TN County Heatmap

Counties are ordered by highest new cases per day (current) based on the most recent 7 day average (per 100k population). Max limit of the color scale set at 200 to increase contrast- some local outbreaks exceed this limit (eg prison outbreaks in Bledsoe, Lake, Wayne Co.)

county_order <-
  tndoh_countynew |> 
  filter(date == max(date)) |> 
  mutate(county_bynewcases = reorder(county, new_cases_per100kpop_7davg)) |> 
  select(county, county_bynewcases)

# no longer interactive due to file size issues.
TN_heatmap_100k <-
  tndoh_countynew |> 
    filter(!county %in% c("Pending", "Out of State")) |>
  left_join(county_order, by = "county") |> 
  mutate(new_cases = replace_na(new_cases, 0)) |> 
  # mutate(lbls =
  #          sprintf("<tr><td><strong><i>%s Co. (%s)</i></strong></tr></td>
  #         <tr><td>Total Cases: %s</tr></td>
  #         <tr><td>New Cases/d/100k, 7dAvg (rank): %1.1f (%s/%s)</tr></td>
  #         <tr><td>New Cases Today (rank): %s (%s/%s)</tr></td>
  #        <tr><td> Test +%% Today: %1.1f%% (%s/%s)</tr></td>
  #         <tr><td>Test +%% (7dAvg): %1.1f%%</tr></td>
  #         <tr><td>Total Deaths: %s</tr></td>", 
  #         toupper(name), paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date)),
  #         comma(total_cases, digits = 0), 
  #         new_cases_per100kpop_7davg, new_cases_per100kpop_7davg_rank, max(new_cases_per100kpop_7davg_rank, na.rm = T), 
  #         comma(new_cases, digits = 0),new_cases_rank, max(new_cases_rank, na.rm = T), 
  #         new_testpos_pct, new_pos_tests, new_tests, 
  #         new_testpos_pct_7davg,
  #         total_deaths) ) |> 
  ggplot(aes(x = date,
             y = county_bynewcases,
             fill = new_cases_per100kpop_7davg)) |> \(x) x +
# TN_heatmap_100k <-
#   TN_heatmap_100k +
  geom_tile(aes(fill = new_cases_per100kpop_7davg,
                color = new_cases_per100kpop_7davg)) +  
  # ggiraph::geom_tile_interactive(aes(tooltip = lbls,
  #                                    data_id = county,
  #                                    fill = new_cases_per100kpop_7davg,
  #                                    color = new_cases_per100kpop_7davg)) +
  viridis::scale_fill_viridis(labels = scales::comma,
                                 limits = c(0, n_ul), oob = scales::squish,
                              option = "plasma") +
  viridis::scale_color_viridis(labels = scales::comma,
                               limits = c(0, n_ul), oob = scales::squish,
                               option = "plasma") +
  labs(fill = "New Cases/100k\n(7dAvg.)",
       color = "New Cases/100k\n(7dAvg.)",
       x = NULL, y = NULL) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b", expand = c(0, 0)) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.background = element_blank(), 
        panel.border = element_blank(),
        legend.title = element_text(size = 14),
        legend.text =  element_text(size = 10))

# x_TN_heatmap_100k <- girafe(ggobj = TN_heatmap_100k,
#             options = list(
#               opts_hover(css = "stroke-width:0.5")))
TN_heatmap_100k
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

Notes:

  • 12/20/2020: Increased upper limit to 200 because 125 just wasn’t enough.

New cases and Hospitalizations by county

Note: assignment of new cases to county usually lags the total count. The State Health Dept. assigns cases to county of residence, rather than testing site. Individual county health departmens may report numbers differently, such as the number of cases diagnosed within the county, regardless of residence

# fig 2 State ####

pl_tnmap_tnsdc <-
  tndoh_countynew |> 
  filter(!is.na(date)) |> 
  filter(date == max(date)) |>
  mutate(total_cases = case_when(total_cases == 0 ~ NA_real_,
                               TRUE ~ total_cases)) |> 
  ggplot() |> \(x) x +
# pl_tnmap_tnsdc <-
#   pl_tnmap_tnsdc +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = as.numeric(total_cases)),
          inherit.aes = F) +
  geom_sf_text(aes(label= comma(total_cases, digits = 0), 
                   geometry = geometry, 
                   color = cut(total_cases, 2)), 
               size = 2.5) + 
  viridis::scale_fill_viridis(labels = scales::comma,
                                 # limits = c(0, 100), oob = scales::squish,
                              option = "plasma") +
    # viridis::scale_fill_viridis(option = "D", direction = -1) +
    scale_color_manual(values = c( "white", "black"), guide = F) +
  labs(fill = "Total Cases:")+
         # caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via TN State Data Center") +
  theme_void() +
  theme(legend.position = "left")  
# ggsave(here::here("figures/pl_tnmap.jpg"), dpi = 600, width = 12, height = 4, units = "in")

pl_tnmap_7dDelta <-
  tndoh_countynew |> #select(new_cases_county_7davg_1wkdelta)
  filter(!is.na(date)) |> 
  filter(date == max(date)) |>
  ggplot() |> \(x) x +
# pl_tnmap_7dDelta <-
#   pl_tnmap_7dDelta +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = new_cases_county_7davg_1wkdelta*100),
          inherit.aes = F) +
  scale_fill_gradient2(low = "green", mid = "white", high = "red") +
  geom_sf_text(aes(label=sprintf("%1.1f",new_cases_county_7davg_1wkdelta*100), geometry = geometry),
               # color = cut(new_cases_county_7davg_1wkdelta, 2)),
               size = 3) +
  # viridis::scale_fill_viridis(option = "D", direction = -1) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "1-Week\nChange %:")+
  # caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via TN State Data Center") +
  theme_void() +
  theme(legend.position = "left")  

pl_tnmap_newhosp <-
  tndoh_countynew |> 
  filter(!is.na(date)) |> 
  filter(date == max(date)) |> 
  ggplot() |> \(x) x +
# pl_tnmap_newhosp <-
#   pl_tnmap_newhosp +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = as.numeric(new_hospitalized)),
          inherit.aes = F) +
  geom_sf_text(aes(label=new_hospitalized, geometry = geometry), 
               color = "black",
               # color = cut(new_hospitalized, 2)), 
               size = 3) + 
  scale_fill_gradient(low = "white", high = "red", 
                      limits = c(0, max(tndoh_countynew$new_hospitalized [tndoh_countynew$date == max(tndoh_countynew$date, na.rm = T)], na.rm = T)))  +
  scale_color_manual(values = c( "white", "black"), guide = F) +
  labs(fill = "New\nHospitalizations:") +
  theme_void() +
  theme(legend.position = "left")  

pl_tnmap_new_tnsdc <-
  tndoh_countynew |> 
  filter(date == max(date)) |> 
  ggplot() |> \(x) x +
# pl_tnmap_new_tnsdc <-
#   pl_tnmap_new_tnsdc +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = as.numeric(new_cases)),
          inherit.aes = F) +
  geom_sf_text(aes(label=new_cases, geometry = geometry), 
               color = "black",
               # color = cut(new_cases, 2)), 
               size = 3) + 
  scale_fill_gradient(low = "white", high = "red", 
                      limits = c(0, max(tndoh_countynew$new_cases[tndoh_countynew$date == max(tndoh_countynew$date, na.rm = T)], na.rm = T)))  +
  labs(fill = "New Cases:") +
       # caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via New York Times") +
  theme_void() +
  theme(legend.position = "left")  
# ggsave(here::here("figures/pl_tnmap_new.jpg"), dpi = 600, width = 12, height = 4, units = "in")


pl_FINAL_tnsdc <-
  cowplot::plot_grid(title_tnsdc,
                   pl_tnmap_new_tnsdc,
                   pl_tnmap_7dDelta,
                   pl_tnmap_newhosp,
                   # pl_tnmap_newtests_tnsdc,
                   pl_tnmap_tnsdc,
                   nrow = 5,
                   rel_heights = c(.1,1,1,1, 1),
                   align = "y")
  unclassified <- 
    tndoh_countynew |> 
    filter(date == max(date),
         county_lower %in% c("pending", "out of state"))

  
    # pl_FINAL_tnsdc
  pl_FINAL_tnsdc +
    annotation_custom(gridExtra::tableGrob(
      tibble(
        Category = unclassified$county,
        Cases = unclassified$new_cases), 
      rows=NULL, cols = NULL,
      theme = gridExtra::ttheme_minimal(base_size = 10)),
      xmin=.85, xmax=.95, ymin=0.66, ymax=.9) +
    annotation_custom(gridExtra::tableGrob(
      tibble(
        Category = unclassified$county,
        Cases = unclassified$new_hospitalized), 
      rows=NULL, cols = NULL,
      theme = gridExtra::ttheme_minimal(base_size = 10)),
      xmin=.85, xmax=.95, ymin=0.31, ymax=.275) +
    annotation_custom(gridExtra::tableGrob(
      tibble(
        Category = unclassified$county,
        Cases = unclassified$total_cases), 
      rows=NULL, cols = NULL,
      theme = gridExtra::ttheme_minimal(base_size = 10)),
      xmin=.85, xmax=.95, ymin=0.01, ymax=.125)
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

Deaths and Crude Case Fatality ratio

The estimate of case fatality is presented as the crude deaths/cases ratio, and is not adjusted for age. This could markedly affect the rate in different counties. This is also affected by the ability to test and identify as many cases as possible, which may differ by county.

pl_deaths <-
  tndoh_countynew |>   
  filter(date == max(date, na.rm = T)) |> 
  ggplot() |> \(x) x +
# pl_deaths <-
#   pl_deaths +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = total_deaths),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.0f", comma(total_deaths, digits = 0)), 
  # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos), 
                   geometry = geometry, 
                   color = cut(total_deaths, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Deaths:") +
  theme_void() +
  theme(legend.position = "left") 


pl_cfr <-
  tndoh_countynew |>   
  filter(date == max(date, na.rm = T)) |> 
  mutate(cfr = total_deaths/total_cases) |> 
  ggplot() |> \(x) x +
# pl_cfr <-
#   pl_cfr +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = cfr*100),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.1f %%", cfr*100), 
  # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos), 
                   geometry = geometry, 
                   color = cut(cfr, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Case Fatality %:") +
  theme_void() +
  theme(legend.position = "left") 




cowplot::plot_grid(title_tnsdc, pl_deaths, pl_tnmap_tnsdc, pl_cfr,
                   nrow = 4, rel_heights = c(.1,1,1, 1))
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 
# ggsave(here::here("figures/pl_deaths_map.jpg"), dpi = 600, width = 12, height = 6, units = "in")

County Population adjusted rates

The estimated Tennessee county populations are taken from the US Census (2019 estimate). The estimate of case fatality is presented as the crude deaths/cases ratio, and is not adjusted for age. This could markedly affect the rate in different counties

pl_deaths <-
  tndoh_countynew |>   
  filter(date == max(date, na.rm = T)) |> 
  ggplot() |> \(x) x +
# pl_deaths <-
#   pl_deaths +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = total_deaths_per100kpop),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.0f", total_deaths_per100kpop), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos), 
                   geometry = geometry, 
                   color = cut(total_deaths_per100kpop, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Deaths per\n100K Pop.:") +
  theme_void() +
  theme(legend.position = "left") 

pl_cases <-
  tndoh_countynew |>   
  filter(date == max(date, na.rm = T)) |> 
  ggplot() |> \(x) x +
# pl_cases <-
#   pl_cases +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = total_cases_per100kpop),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.1f%%", total_cases_per100kpop/1000), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos), 
                   geometry = geometry, 
                   color = cut(total_cases_per100kpop, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Population % Affected\n(Cases %Pop.):") +
  theme_void() +
  theme(legend.position = "left") 

pl_newcases <-
  tndoh_countynew |> #select(new_cases_per100kpop_7davg)   
  filter(date == max(date, na.rm = T)) |> 
  ggplot() |> \(x) x +
# pl_newcases <-
#   pl_newcases +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = new_cases_per100kpop_7davg),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.0f", new_cases_per100kpop_7davg), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos), 
                   geometry = geometry, 
                   color = cut(new_cases_per100kpop_7davg, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Daily Cases per\n100K Pop.\n(7d Avg):") +
  theme_void() +
  theme(legend.position = "left") 

pl_cfr <-
  tndoh_countynew |> 
    mutate(total_tests_pctpop = total_tests/population*100000) |> 
    # select(total_tests_per100kpop)   
  filter(date == max(date, na.rm = T)) |> 
  ggplot() |> \(x) x +
# pl_cfr <-
#   pl_cfr +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = total_tests_pctpop),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.1f ", total_tests_pctpop/100), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos), 
                   geometry = geometry, 
                   color = cut(total_tests_pctpop, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Tests per\n100k Population:") +
  theme_void() +
  theme(legend.position = "left") 

cowplot::plot_grid(  title_tnsdc, pl_newcases,pl_deaths, pl_cases,pl_cfr,
                   nrow = 5, rel_heights = c(.1,1,1,1, 1))
<<<<<<< HEAD

=======

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e 
# ggsave(here::here("figures/pl_deaths_map.jpg"), dpi = 600, width = 12, height = 6, units = "in")

US Data from COVID Tracking Project

Testing and Case data for these maps are taken from the COVID Tracking project, and the 2019 state population estimates are taken from the United States Census Bureau.

New Cases

mouse over state for details

max_newcaserate <- max(covus_currentdata_2163$new_cases_state_7davg_100k[covus_currentdata_2163$date == max(covus_currentdata_2163$date)], na.rm = T)
                                                  
pl_us_cases_100k <-
  covus_currentdata_2163 |> 
  filter(date == max(date)) |> 
  ggplot(aes(geometry = geometry,
             fill = new_cases_state_7davg_100k)) |> \(x) x +
# pl_us_cases_100k <-
#   pl_us_cases_100k +
  ggiraph::geom_sf_interactive(aes(tooltip = tt,
                                   data_id = state), size = .2) +
  viridis::scale_fill_viridis(labels = scales::comma,
                               limits = c(0, max_newcaserate), oob = scales::squish,
                              option = "plasma") +
  labs(fill = "Cases/100k/d\n(7dAvg)") +
  theme_void() 



ggiraph::girafe(code = print( pl_us_cases_100k))
<<<<<<< HEAD
=======
>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

US Heatmap

max_newcaserate_alldates <- max(covus_currentdata_2163$new_cases_state_7davg_100k, na.rm = T)

state_order <- 
  covus_currentdata_2163 |> 
  filter(date == max(date)) |> 
  mutate(state_bynewcases = reorder(state, new_cases_state_7davg_100k)) |> 
  select(state, state_bynewcases)

us_heatmap_100k <-
  covus_currentdata_2163 |> 
  left_join(state_order, by = "state") |> 
  filter(date >= lubridate::ymd("2020-04-01")) |> 
  ggplot(aes(x = date,
             y = state_bynewcases)) |> \(x) x +
# us_heatmap_100k <-
#   us_heatmap_100k +
  geom_tile(aes(fill = new_cases_state_7davg_100k,
                color = new_cases_state_7davg_100k)) +
  # ggiraph::geom_tile_interactive(aes(tooltip = tt,
  #                                    data_id = state,
  #                                    fill = new_cases_state_7davg_100k,
  #                                    color = new_cases_state_7davg_100k)) +
  viridis::scale_fill_viridis(labels = scales::comma,
                               limits = c(0, max_newcaserate_alldates), oob = scales::squish,
                              option = "plasma") +
    viridis::scale_color_viridis(labels = scales::comma,
                               limits = c(0, max_newcaserate_alldates), oob = scales::squish,
                              option = "plasma") +
  labs(fill = "Cases/100k/d\n(7dAvg)",
       color  = "Cases/100k/d\n(7dAvg)",
       x = NULL, y = NULL) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b", expand = c(0, 0)) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 12),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.background = element_blank(), 
        panel.border = element_blank(),
        legend.title = element_text(size = 14),
        legend.text =  element_text(size = 12))
us_heatmap_100k

You made it to the end, so here’s a Dad joke:

<<<<<<< HEAD

I used to hate facial hair…but then it grew on me.

=======

What do an apple and an orange have in common?

Neither one can drive.

>>>>>>> 89a401a8efd7c49dd29e5f9169754337f657961e

Session Info

Hmisc::markupSpecs$html$session()
 R Under development (unstable) (2020-12-09 r79601)
 Platform: x86_64-w64-mingw32/x64 (64-bit)
 Running under: Windows 10 x64 (build 17763)
 
 Matrix products: default
 
 attached base packages:
 [1] stats     graphics  grDevices utils     datasets  methods   base     
 
 other attached packages:
  [1] gdtools_0.2.2       ggiraph_0.7.8       formattable_0.2.0.1
  [4] leafpop_0.0.6       leaflet_2.0.3       sf_0.9-6           
  [7] forcats_0.5.0       stringr_1.4.0       dplyr_1.0.2        
 [10] purrr_0.3.4         readr_1.4.0         tidyr_1.1.2        
 [13] tibble_3.0.4        ggplot2_3.3.2       tidyverse_1.3.0
To cite R in publication use:

R Core Team (2020). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.